# ----------------------------------------------------------------------
# How Does Shaming Human Rights Violators Abroad Shape Attitudes at Home?
# Supplementary Analyses (observational data)
# ----------------------------------------------------------------------

# ----------------------------------------------------------------------
# load all relevant themes and packages
# ----------------------------------------------------------------------
rm(list=ls())
library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("RItools")
library("estimatr")
library("fastDummies")
library("ggpubr")
library("car")
library("table1")#
library("psy")#
library("dplyr")
library("ggstatsplot")
library("stringr")
library("haven")
library("lfe")
library("naniar")
library("fixest")
library("Hmisc")
library("quanteda")
library("readtext")
library("htmltools")
library("devtools")
library("sjPlot")
library("kableExtra")
library("maps")
library("countrycode")
library("multiwayvcov")
library("lmtest")
library("texreg")

# ----------------------------------------------------------------------
# Create Figure 1 (Heatmap)
# ----------------------------------------------------------------------

#download `dataset1.csv` (UPR shaming dataset)
merged_shaming <- read_csv("dataset1.csv")

#create variable of average of shaming per country
merged_shaming<-
  merged_shaming %>%
  group_by(country) %>%
  mutate(.,
         mean_recs = mean(all_recs))

vars <- c("country_year","country","mean_recs")
maps <- merged_shaming[vars]
maps$country_year<-NULL
maps <-maps %>% distinct()

# load map
world_map <- map_data("world")
world_map$country<-countrycode(sourcevar = world_map$region, origin = "country.name", 
                               destination = "cowc")

count_map <- left_join(world_map,maps, by = c("country"))
count_map$mean_recs<-ifelse(is.na(count_map$mean_recs),0,count_map$mean_recs)

count_map <-
  count_map %>% 
  mutate(.,
         Shaming = case_when(
           mean_recs > 150 ~ "Over 150 average \n UPR recommendations",
           mean_recs > 100 & mean_recs < 150 ~ "100-150 average \n UPR recommendations",
           mean_recs > 50 & mean_recs < 100 ~ "50-100 average \n UPR recommendations",
           mean_recs > 10 & mean_recs < 50 ~ "10-50 average \n UPR recommendations",
           mean_recs < 10 ~ "Under 10 average \n UPR recommendations"
         ))

ggplot(count_map, aes(long, lat, group = group))+
  geom_polygon(aes(fill = Shaming), color = "white", 
               size =.01)+
  labs(title = "Global Distribution of Shamers in the UPR \n (2008-2020)",
       fill = "",
       caption = "")+
  scale_fill_grey(limits = c( "Over 150 average \n UPR recommendations",
                              "100-150 average \n UPR recommendations",
                              "50-100 average \n UPR recommendations",
                              "10-50 average \n UPR recommendations",
                              "Under 10 average \n UPR recommendations"), na.value = "gray", na.translate = F) +
  theme_void(base_size = 10)+
  theme(legend.position="bottom",
        text = element_text(family = "serif"),
        plot.title = element_text(hjust = 0.5, face = "bold"))

ggsave("../../Main_Figs/fig1.pdf", width = 8, height = 6)

#who are the countries that shame the most? (on average)
maps%>%
  group_by(country) %>% 
  tally(mean_recs) %>% 
  top_n(n = 5)

#France (190), Spain (173), Canada (158), Mexico (133) and Norway (125)

# ----------------------------------------------------------------------
# Observational Analyses
# ----------------------------------------------------------------------

#load `dataset2.csv` (CSES and UPR shaming)
dataset2 <- read_csv("dataset2.csv")

#Model_1a: Main model OLS
#Independent variable -- shaming UPR (harsh categories) per country per year (t-1)
#Dependent variable -- Incumbent vote
#Fixed effects -- Country and year
#Errors clustered at the country-year level

Model_1a<- felm(incumbent_vote ~ lag_recs
                |year + country.x | 0 |country_year,
                data = dataset2)

summary(Model_1a)

#Model_2a: Main model with individual-level control
#governments who care about HRs are more likely to shame, and respondents who care about HRs are also more likely to vote for them, regardless of shaming. 

Model_2a <- felm(incumbent_vote ~ lag_recs + imputed_iddistance +
                   imputed_age + as.factor(imputed_gender) +
                   imputed_education 
                 |country.x + year | 0 |country_year,
                 data = dataset2)

summary(Model_2a) 

#Model_3a: Main model with physical integrity rights, state capacity, and demographic controls
Model_3a <- felm(incumbent_vote ~ lag_recs
                 + imputed_iddistance
                 + imputed_theta 
                 + imputed_capacity +
                   imputed_age + as.factor(imputed_gender) +
                   imputed_education
                 |country.x + year | 0 |country_year,
                 data = dataset2)

summary(Model_3a)

# Create table
tab1<-stargazer(Model_1a,Model_2a,Model_3a,
                out= "../Tables/SA/C16.tex",  style="qje",
                title="The effect of UPR Shaming on incumbent voting (CSES)",
                keep = c("lag_recs"),
                dep.var.labels = "Incumbent vote",
                covariate.labels = "UPR Shaming (t-1)",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.1, 0.05, 0.01),
                label = "tab:H1_OLS",
                #star.char = c("", "", ""),
                add.lines = list(
                  c("Individual-level Controls",  "No", "Yes", "Yes"),
                  c("Physical Integrity (t-1)",  "No", "No", "Yes"),
                  c("State Capacity (t-1)",  "No", "No", "Yes"),
                  c("Year FEs", "Yes", "Yes", "Yes"),
                  c("Country FEs",  "Yes", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Robust standard errors are clustered at the country-year level. Missing data for control variables is imputed (average value by country-year)."))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Robust standard errors are clustered
at the country-year level. Missing data for control variables is imputed 
(average value by country-year).
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Tables/SA/C16.tex')

### OLS models - alt. clustering (country) ####

#Model_1A: OLS country-level clustering
#Country, and year fixed effects
#Errors clustered at the country level
Model_1A<- felm(incumbent_vote ~ lag_recs
                |year + country.x | 0 |country.x,
                data = dataset2)

summary(Model_1A)

#Model_2A: Model_1A with demographic controls

Model_2A<- felm(incumbent_vote ~ lag_recs + imputed_iddistance +
                  imputed_age + as.factor(imputed_gender) +
                  imputed_education 
                |country.x + year  | 0 |country.x,
                data = dataset2)

summary(Model_2A) 

#Model_3A: Model_1A with physical integrity rights and demographic controls
Model_3A <- felm(incumbent_vote ~ lag_recs
                 + imputed_theta
                 + imputed_capacity
                 + imputed_iddistance +
                   imputed_age + as.factor(imputed_gender) +
                   imputed_education
                 |country.x + year | 0 |country.x,
                 data = dataset2)

summary(Model_3A)

#Alternative modeling of independent variable (ALL recommendations)
Model_A<- felm(incumbent_vote ~ lag_all
               |year + country.x | 0 |country_year,
               data = dataset2)

summary(Model_A)

#Model_B: Main model with individual-level control
Model_B <- felm(incumbent_vote ~ lag_all + imputed_iddistance +
                  imputed_age + as.factor(imputed_gender) +
                  imputed_education 
                |country.x + year | 0 |country_year,
                data = dataset2)

summary(Model_B) 

#Model_C: Main model with physical integrity rights and demographic controls
Model_C <- felm(incumbent_vote ~ lag_all
                + imputed_iddistance
                + imputed_theta +
                  imputed_capacity +
                  imputed_age + as.factor(imputed_gender) +
                  imputed_education
                |country.x + year | 0 |country_year,
                data = dataset2)

summary(Model_C)


# Create table
tab1<-stargazer(Model_2A,Model_2A,Model_3A,
                Model_A,Model_B,Model_C,
                out= "../Tables/SA/C20.tex",  style="qje",
                title="Alternative Modeling (CSES)",
                keep = c("lag_recs","lag_all"),
                dep.var.labels = "Incumbent vote",
                covariate.labels = c("UPR Shaming (t-1)","ALL UPR Recommendations (t-1)"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.1, 0.05, 0.01),
                label = "tab:H1_alt1",
                # star.char = c("", "", ""),
                add.lines = list(
                  c("Individual-level Controls",  "No", "Yes", "Yes", "No", "Yes", "Yes"),
                  c("Physical Integrity (t-1)",  "No", "No", "Yes","No", "No", "Yes"),
                  c("State Capacity (t-1)",  "No", "No", "Yes","No", "No", "Yes"),
                  c("Year FEs", "Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                  c("Country FEs",  "Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                  c("Clustering at country level", "Yes", "Yes", "Yes", "No", "No", "No")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE,
                notes.align = "l",
                notes = c("Missing data for control variables is imputed (average value by country-year)."))

note.latex <- "\\multicolumn{7}{c} {\\parbox[t]{12cm}{ \\textit{Notes:}Missing data for 
                control variables is imputed (average value by country-year). 
                * $p < .1$, ** $p < .05$, *** $p < .01$ }} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Tables/SA/C20.tex')


### GLM models ####

#Model_1b: GLM models
#Country, and year fixed effects
#Errors clustered at the country-year level

Model_1b<- feglm(incumbent_vote ~ lag_recs
                 |year + country.x, cluster = "country_year",
                 data = dataset2)

summary(Model_1b)

#Model_2b: Model_1b with demographic controls
Model_2b<- feglm(incumbent_vote ~ lag_recs + imputed_iddistance +
                   imputed_age + as.factor(imputed_gender) +
                   imputed_education 
                 |year + country.x, cluster = "country_year",
                 data = dataset2)

summary(Model_2b) 


#Model_3b: Model_1b with physical integrity rights and demographic controls
Model_3b <- feglm(incumbent_vote ~ lag_recs
                  + imputed_theta
                  + imputed_capacity
                  + imputed_iddistance +
                    imputed_age + as.factor(imputed_gender) +
                    imputed_education
                  |year + country.x, cluster = "country_year",
                  data = dataset2)

summary(Model_3b)

# Create table
texreg(list(Model_1b,Model_2b,Model_3b), 
       dcolumn = TRUE, booktabs = TRUE, use.packages = FALSE, 
       label = "tab:H1_GLM", 
       caption = "The effect of UPR shaming on incumbent voting (CSES). Logistic Regression Models. Robust standard errors are clustered at the country-year level. Missing data for control variables is imputed (average value by country-year).", 
       stars = c(0.1,0.05,0.001), digits = 4, 
       float.pos = "h",
       custom.coef.map = list("lag_recs" = "UPR Shaming (t-1)"),
       custom.gof.rows = list("Individual-level Controls" = c("No", "Yes", "Yes"),
                              "Physical Integrity (t-1)" = c("No", "No", "Yes"),
                              "Imputed Capacity (t-1)" = c("No", "No", "Yes"),
                              "Year FEs" = c("Yes", "Yes", "Yes"),
                              "Country FEs" = c("Yes", "Yes", "Yes")),
       file = "../Tables/SA/C19.tex")


###########################
# Analysis -- H2 (WVS)
###########################

#load `dataset3.csv` (WVS and UPR shaming)
dataset3 <- read_csv("dataset3.csv")

### Main OLS models ####

#Model_1c: Main model OLS
#Independent variable -- "harsh" recommendations
#Country, and year fixed effects
#Errors clustered at the country-year level

Model_1c<-felm(HR_perception ~ lag_recs
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_1c)

##Model_2c: Add individual-level controls
Model_2c<-felm(HR_perception ~ lag_recs  + as.factor(imputed_sex)
               + imputed_age + imputed_education + imputed_id  + imputed_income
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_2c)


##Model_3c: Add physical integrity + state capacity 
Model_3c<-felm(HR_perception ~ lag_recs  + as.factor(imputed_sex) 
               + imputed_age + imputed_education + imputed_id + imputed_income + imputed_theta
               + imputed_capacity
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_3c)


# Create table
tab1<-stargazer(Model_1c,Model_2c,Model_3c,
                out= "../Tables/SA/C17.tex",  style="qje",
                title="The effect of UPR shaing on perceptions of human rights (WVS)",
                keep = "lag_recs",
                dep.var.labels = "Human Rights Perceptions",
                covariate.labels = "UPR Shaming (t-1)",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.1,0.05,0.01),
                label = "tab:H3_OLS",
                #star.char = c("", "", ""),
                add.lines = list(
                  c("Individual-level Controls",  "No", "Yes", "Yes"),
                  c("Physical Integrity (t-1)",  "No", "No", "Yes"),
                  c("State Capacity (t-1)",  "No", "No", "Yes"),
                  c("Year FEs", "Yes", "Yes", "Yes"),
                  c("Country FEs",  "Yes", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Robust standard errors are clustered at the country-year level. Missing data for control variables is imputed (average value by country-year)."))

note.latex <- "\\multicolumn{4}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Robust standard errors 
are clustered at the country-year level. Missing data for control variables is imputed (average 
value by country-year). 
* $p < .1$, ** $p < .05$, *** $p < .01$ }} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Tables/SA/C17.tex')


### Alt. clustering (country) ####

#Model_1C: OLS country-level clustering
#Country, and year fixed effects
#Errors clustered at the country level

Model_1C<-felm(HR_perception ~ lag_recs
               |country + year| 0 |country,
               data = dataset3)

summary(Model_1C)


##Model_2C: Add demographic controls
Model_2C<-felm(HR_perception ~ lag_recs  + as.factor(imputed_sex)
               + imputed_age + imputed_education + imputed_id + imputed_income
               |country + year| 0 |country,
               data = dataset3)

summary(Model_2C)

##Model_3C: Add physical integrity + capacity
Model_3C<-felm(HR_perception ~ lag_recs  + as.factor(imputed_sex) + imputed_theta + imputed_capacity
               + imputed_age + imputed_education + imputed_id + imputed_income
               |country + year| 0 |country,
               data = dataset3)

summary(Model_3C)


#Alternative modeling of independent variable as ALL recommendations
Model_1D<-felm(HR_perception ~ lag_all
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_1D)


##Model_2C_all: Add demographic controls
Model_2D<-felm(HR_perception ~ lag_all  + as.factor(imputed_sex)
               + imputed_age + imputed_education + imputed_id + imputed_income
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_2D)

##Model_3C_all: Add physical integrity + state capacity
Model_3D<-felm(HR_perception ~ lag_all  + as.factor(imputed_sex) + imputed_theta + imputed_capacity
               + imputed_age + imputed_education + imputed_id + imputed_income
               |country + year| 0 |country_year,
               data = dataset3)

summary(Model_3D)

# Create table
tab1<-stargazer(Model_1C,Model_2C,Model_3C,
                Model_1D, Model_2D, Model_3D,
                out= "../Tables/SA/C21.tex",  style="qje",
                title="Alternative modeling (WVS)",
                keep = c("lag_recs", "lag_all"),
                dep.var.labels = "Human Rights Perceptions",
                covariate.labels = c("UPR Shaming (t-1)", "All UPR Recommendations (t-1)"), 
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.11,0.05,0.01),
                label = "tab:H3_alt1",
                # star.char = c("", "", ""),
                add.lines = list(
                  c("Individual-level Controls",  "No", "Yes", "Yes","No", "Yes", "Yes"),
                  c("Physical Integrity (t-1)",  "No", "No", "Yes","No", "No", "Yes"),
                  c("State Capacity (t-1)",  "No", "No", "Yes","No", "No", "Yes"),
                  c("Year FEs", "Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                  c("Country FEs",  "Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                  c("Clustering at country level", "Yes", "Yes", "Yes", "No", "No", "No")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE,
                notes.align = "l",
                notes = c("Missing data for control variables is imputed (average value by country-year)."))

note.latex <- "\\multicolumn{7}{c} {\\parbox[t]{12cm}{ \\textit{Notes:}Missing data 
for control variables is imputed (average value by country-year).}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Tables/SA/C21.tex')

###########################
# Placebo Tests
###########################

#We regress our IV (shaming) at time t+1 over our DV at time t to reduce concerns of selection bias

#H1 Model

Placebo_H1<- felm(incumbent_vote ~ lead_recs
                  |year + country.x | 0 |country_year,
                  data = dataset2)

summary(Placebo_H1)

#H2 Model

Placebo_H2<-felm(HR_perception ~ lead_recs
                 |country + year| 0 |country_year,
                 data = dataset3)

summary(Placebo_H2)


# Create table
tab1<-stargazer(Placebo_H1,Placebo_H2,
                out= "../Tables/SA/C18.tex",  style="qje",
                title="Placebo test: the effect of UPR shaming at t+1 on dependent variables at t",
                keep = c("lead_recs"),
                dep.var.labels = c("Incumbent vote","Human rights perceptions"),
                covariate.labels = "UPR Shaming (t+1)",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = c(0.1,0.05,0.01),
                label = "tab:Placebos",
                #star.char = c("", "", ""),
                add.lines = list(
                  c("Year FEs", "Yes", "Yes"),
                  c("Country FEs",  "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c("Robust standard errors are clustered at the country-year level."))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Robust standard errors are clustered at the country-year level.}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Tables/SA/C18.tex')

###########################
# Descriptive Statistics
###########################

#Create descriptive statistics tables

#H1 descriptive
dataset2 <- dataset2 %>% 
  mutate(.,
         female = case_when(imputed_gender==1~1,
                            imputed_gender==2~0))

disc_table <- select(dataset2,
                     c(incumbent_vote, lag_recs, imputed_theta, imputed_iddistance, 
                       imputed_age, female, imputed_education))

stargazer(as.data.frame(disc_table),
          covariate.labels = c("Incumbent Vote","UPR Shaming", "Physical integrity rights", "Ideological distance",
                               "Age", "Gender (Female)", "Education"),
          title = "Descriptive Statistics - CSES (H1)",
          label = "tab:dsc_H1",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="../Tables/SA/C22.tex")


#H2 descriptives
dataset3 <- dataset3 %>% 
  mutate(.,
         female = case_when(imputed_sex==1~1,
                            imputed_sex==2~0))

disc_table <- select(dataset3,
                     c(HR_perception, lag_recs, imputed_theta, 
                       imputed_id, imputed_age, female, imputed_education, imputed_income))

stargazer(as.data.frame(disc_table),
          covariate.labels = c("Human Rights Perceptions","UPR Shaming", "Physical integrity rights", "Ideology",
                               "Age", "Gender (Female)", "Education", "Income"),
          title = "Descriptive Statistics - WVS (H2)",
          label = "tab:dsc_H3",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="../Tables/SA/C23.tex")

# ----------------------------------------------------------------------
# Shaming awareness analysis
# ----------------------------------------------------------------------

### read dataset4
data <- read_csv("dataset4.csv")
data <- data[3:nrow(data),]

# remove inattentive
data_attentive <- data %>% 
  filter(., attention==2)

#how many heard of criticism?

(nrow(data_attentive %>% 
  filter(., HR_month==1))) / nrow(data_attentive) 

#86.3% say they heard in last month

#Create table C15

data_attentive$age<-as.numeric(data_attentive$age)
data_attentive$education<-as.numeric(data_attentive$education)
  
data_attentive <- data_attentive %>%
  mutate(.,
         HR_month2 = case_when(
           HR_month == "1" ~ 1,
           HR_month == "2" ~ 0),
         Female = case_when(
           gender == "1" ~ 1,
           gender == "2" ~ 0),
         Male = case_when(
           gender == "1" ~ 0,
           gender == "2" ~ 1),
         Democrat = case_when(
           partisan == "1" ~ 1,
           partisan == "2" ~ 1,
           partisan == "3" ~ 1,
           partisan == "4" ~ 0,
           partisan == "5" ~ 0,
           partisan == "6" ~ 0,
           partisan == "7" ~ 0),
         Independent = case_when(
           partisan == "1" ~ 0,
           partisan == "2" ~ 0,
           partisan == "3" ~ 0,
           partisan == "4" ~ 1,
           partisan == "5" ~ 0,
           partisan == "6" ~ 0,
           partisan == "7" ~ 0),
         Republican = case_when(
           partisan == "1" ~ 0,
           partisan == "2" ~ 0,
           partisan == "3" ~ 0,
           partisan == "4" ~ 0,
           partisan == "5" ~ 1,
           partisan == "6" ~ 1,
           partisan == "7" ~ 1))

disc_table <- dplyr::select(data_attentive, c(HR_month2, Female,
                                              Male, Democrat, Independent,
                                              Republican,age, education))

#add in more demographics
stargazer(as.data.frame(disc_table),
          covariate.labels = c("Recall HRs shaming", "Female", "Male", 
                               "Democrat", "Independent", "Republican", "Age", "Education"),
          
          title = "Descriptive Statistics",
          label = "tab:dsc_follow",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="../Tables/SA/C15.tex")


